#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <assert.h>

#include "myrandom.h"
#include "func.h"

static double h(double x,double alpha);


/* Generacion de una variable Exponencial de parametro l. */
/* Metodo de la trasformada inversa. */
double E(double l){
	double u = 0.0;
	assert(l != 0);
	u = (1.0/(float)ULONG_MAX) * mzran13();
	return  (-1.0/l)*log(u);
}


/* Generacion de una variable Geometrica de parametro p. */
/* Metodo de la trasformada inversa. */
unsigned int G(double p) {
	double U = 0.0, res = 0.0;
      
	U = (1.0/(float)ULONG_MAX) * mzran13();
	res = (int)(log(U)/log(1.0-p)) + 1;
  
	return res;
}

/*Generacion de v.a Normal con el metodo polar*/
double *Normal_MP(double mu, double sigma){
	double *res = NULL, u1 = 0.0, u2 = 0.0, v1 = 0.0, v2 = 0.0, s = 2.0;
	
	res = (double*)calloc(2, sizeof(double));
	
	while (s > 1.0){
		u1 = (1.0/(float)ULONG_MAX) * mzran13();
		u2 = (1.0/(float)ULONG_MAX) * mzran13();
		v1 = (u1*2) - 1.0;
		v2 = (u2*2) - 1.0; 
		s = pow(v1, 2) + pow(v2, 2);
	}
	
	res[0] = ((pow((((-1)*(2)*log(s))/s), 0.5)*v1)*sigma)+mu;
	res[1] = ((pow((((-1)*(2)*log(s))/s), 0.5)*v2)*sigma)+mu;
	
	return res;
}

/*Generacion de una lognormal con una normal estandar*/
double lognormal(double mu, double sigma){
	double *normal = Normal_MP(0.0, 1.0);
	double res = normal[0];
	free(normal);
	
	return (exp(mu + sigma * res));
}

/*Generacion de una gamma*/
/*Metodo de aceptacion y rechazo. Rechaza con una exponencial de igual media*/
double gamma(double alpha, double beta){
	double alpha_aux = 1.0/alpha;
	double c = h((alpha-1.0)/(-1.0*(alpha_aux-1.0)), alpha);
	double u=0.0, v = 0.0, y = 0.0;
	
	do{
		u = (1.0/(float)ULONG_MAX) * mzran13();
		v = (1.0/(float)ULONG_MAX) * mzran13();
	
		y = -alpha*log(v);
		
	}while (u>= h(y,alpha)/c);
	
	return beta*y;
}

static double h(double x,double alpha){
	
	double alpha_aux = 1.0/alpha;
	double k = 1.0/(exp(gammln(alpha)));
	
	return (k/alpha_aux) * exp(x*(alpha_aux - 1.0)) * pow(x, alpha-1.0);
}
	
	

